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ABSTRACT 

We design an algorithm for computing the p-curvature of a 
differential system in positive characteristic p. For a system 
of dimension r with coefficients of degree at most d, its com¬ 
plexity is 0~{pdr‘^) operations in the ground field (where 
ui denotes the exponent of matrix multiplication), whereas 
the size of the output is about pdr^. Our algorithm is then 
quasi-optimal assuming that matrix multiplication is {i.e. 
LO = 2). The main theoretical input we are using is the exis¬ 
tence of a well-suited ring of series with divided powers for 
which an analogue of the Cauchy-Lipschitz Theorem holds. 
Categories and Subject Descriptors: 

1.1.2 [Computing Methodologies]: Symbolic and Alge¬ 
braic Manipulation - Algebraic Algorithms 

Keywords: Algorithms, complexity, differential equations, 
p-curvature. 

1. INTRODUCTION 

We study in this article algorithmic questions related to 
linear differential systems in positive characteristic. Let k 
be an arbitrary field of prime characteristic p, and A be 
an r X r matrix with entries in the field k{x) of rational 
functions over k. A simple-to-define, yet very important 
object attached to the differential system Y' = AV is its 
so-called p-curvature. It is the p-th iterate d\ of the map 
Oa ■ k{xY —>■ k{xY that sends v to v' — Av. It turns out 
that it is fc(a;)-linear. It is moreover classical that its matrix 
with respect to the canonical basis of k{xY is equal to the 
term Ap of the recursive sequence {Ai)i defined by 

Ai = —A and Ai+i = A'^ — A ■ Ai for i > 1. (1) 

In all what follows, we will thus deliberately identify the 
matrix Ap with the p-curvature of Y' — AY. The above 
recurrence yields an algorithm for computing it, sometimes 
referred to as Katz’s algorithm. 

The p-curvature is related to solutions; it measures to 
what extent the usual Cauchy-Lipschitz theorem applies in 
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characteristic p. More precisely, at an ordinary point, the 
system Y' = AY admits a fundamental matrix of power 
series solutions in fc[[a;]] if and only if the p-curvature Ap 
vanishes. In this case, the system Y' = AY even admits a 
fundamental matrix of solutions which are rational functions 
in k(x). More generally, the dimension of the kernel of Ap 
is equal to the dimension of the space of rational function 
solutions of Y' = AY. 

The primary importance of the notion of p-curvature re¬ 
lies in its occurrence in one of the versions of the celebrated 
Grothendieck-Katz conjecture [19, 20, 12, 30). This con¬ 
jecture, first formulated by Alexandre Grothendieck in the 
late 1960s, is a local-global principle for linear differential 
systems, which states that a linear differential system with 
rational function coefficients over a function field admits a 
fundamental matrix of algebraic solutions if and only if its 
p-curvatures vanish for almost all primes p. 

In computer algebra, p-curvature has been introduced by 
van der Put [22, 23], who popularized it as a tool for fac¬ 
toring differential operators in characteristic p. Gluzeau [13] 
generalized the approach to the decomposition of differen¬ 
tial systems over k{x). The p-curvature has also been used 
by Gluzeau and van Hoeij [14] as an algorithmic filter for 
computing exponential solutions of differential operators in 
characteristic zero. 

Gomputing efficiently the p-curvature is in itself a chal¬ 
lenging problem, especially for large values of p. Our initial 
motivation for studying this question emerged from concrete 
applications, in lattice path combinatorics ]6, 7] and in sta¬ 
tistical physics [3]. In this article, we address the question of 
the computation of Ap in good complexity, with respect to 
three parameters: the dimension r of the system Y' = AY, 
the maximum degree d of the rational function entries of 
A, and the characteristic p of the ground field. In terms of 
these quantities, the arithmetical size of Ap is generically 
proportional to pdr^ if r > 1. 

Previous work. Gluzeau [13, Prop. 3.2] observed that 
the direct algorithm based on recurrence (1) has complexity 
0~(p^dr“), where oj is the matrix multiplication exponent 
and the soft-O notation 0~{) hides polylogarithmic factors. 
Gompared to the size of the p-curvature, this cost is good 
with respect to r and d, but not to p. The first subquadratic 
algorithm in p, of complexity 0~(p^''"“^^), was designed in [9, 
§6.3]. In some special cases, additional partial results were 
obtained in [9], notably an algorithm of quasi-linear cost 
0~(p) for certain systems of order r = 2. However, the ques¬ 
tion of designing a general algorithm for computing Ap with 
quasi-linear complexity in p remained open. In a related, but 


different direction, the article [4] proposed an algorithm for 
compnting the characteristic polynomial of the p-curvature 
in time essentially linear in y/p, without computing Ap itself. 

Contribution. We prove that the p-curvature Ap can be 
computed in quasi-linear time with respect to p. More pre¬ 
cisely, our main result (Theorem 4.2) states that 0''(pdr“) 
operations in k are sufficient for this task. This complexity 
result is quasi-optimal not only with respect to the main pa¬ 
rameter p, but also to d; with respect to the dimension r, it 
is as optimal as matrix multiplication. Moreover the algo¬ 
rithm we obtain is highly parallelizable by design. The key 
tools underlying the proof of Theorem 4.2 are the notion of 
divided power rings in characteristic p, and a new formula 
for the p-curvature (Propositions 4.3 and 4.4) in terms of 
divided power series. Crucial ingredients are the fact that 
a Cauchy-Lipschitz theorem for differential systems holds 
over divided power rings (Proposition 3.4) and the fact that 
Newton iteration can be used to efhciently compute (trunca¬ 
tions of) fundamental matrices of divided power solutions. 

Structure of the paper. In Section 2, we recall the main 
theoretical properties of the basic objects used in this article. 
Section 3 is devoted to the existence and the computation 
of solutions of differential systems in divided power rings. 
In Section 4, we move to the main objective of the arti¬ 
cle, the computation of the p-curvature: after relating Ap 
to the framework of divided powers, we describe our main 
algorithm for Ap, of complexity 0~{pdr‘^). We conclude in 
Section 5 by describing the implementation of our algorithm 
and some benchmarks. 

Complexity measures. Throughout this article, we es¬ 
timate the cost of our algorithms by counting arithmetic 
operations in the base ring or field at unit cost. 

We use standard complexity notations. The letter lu refers 
to a feasible exponent for matrix multiplication {i.e. there 
exists an algorithm for multiplying n x n matrices over a 
ring 21 with at most 0{n'^) operations in 21); the best known 
bound is a; < 2.3729 from [15]. The soft-O notation O'(-) 
indicates that poly logarithmic factors are omitted; in par¬ 
ticular, we will use the fact that many arithmetic operations 
on univariate polynomials of degree d can be done in 0~{d) 
operations: addition, multiplication, Chinese remaindering, 
etc, the key to these results being fast polynomial multiplica¬ 
tion [27, 26, 11, 18]. A general reference for these questions 
is [17]. 

2. THEORETICAL SETTING 

We introduce and briefly recall the main properties of the 
theoretical objects we are going to use in this article. All 
the material presented in this section is classical; a general 
reference is [24]. 

Definitions and notations. Let 21 be a commutative ring 
with unit. We recall that a derivation on 21 is an additive 
map ' : 21 — >■ 21, satisfying the Leibniz rule (fg)' = f'g + fg' 
for all f,g£Ql. The image /' of / under the derivation is 
called the derivative of /. From now on, we assume that 
21 is equipped with a derivation. A differential system with 
coefficients in 21 is an equation of the form Y' = AY where 
A is a given r xr matrix with coefficients in 21 (for a certain 
positive integer r), the unknown T is a column vector of 
length r and Y' denotes the vector obtained from Y by tak¬ 
ing the derivative component-wise. The integer r is called 


the dimension of the system. We recall briefly that a linear 
differential equation: 

-I- ■ • • -I- aiy' + aoy = 0 (with Oi € 21) (2) 


can be viewed as a particular case of a differential system. 
Indeed, defining the companion matrix 


C = 


/ 

1 



V 1 -V/ 


(3) 


and A = *(7, the solutions of the system Y' = AY are 
exactly the vectors of the form ^{y,y', ■ ■ ■ where y 

is a solution of (2). In this correspondence, the order of 
the differential equation agrees with the dimension of the 
associated differential system. 


Differential modules. A differential module over 21 is a 
pair {M, d) where M is an 21-module and d : M —>■ M is an 
additive map satisfying a Leibniz-like rule, which is: 

yfe‘2l,yxeM, d{fx) = f ■ x +f ■ d{x). (4) 


There exists a canonical one-to-one correspondence between 
differential systems and differential modules (M, d) for which 
M = 21’’ for some r: to a differential system Y' = AV of di¬ 
mension r, we attach the differential module (21’", Oa) where 
Sa : —>■ is the function mapping X to X' — AX. Under 
this correspondence, the solutions of Y' = AY are exactly 
vectors in the kernel of Oa- 

To a differential equation as (2), one can associate the 
differential operator L — Ord’’ -I- a, — id’’”^ -I- • • • -I- aid -I- ao; 
it lies in the non-commutative ring 21(d), endowed with the 
usual addition of polynomials and a multiplication ruled by 
the relation d- f = f -8 + f for all / G 21 (note that, as often 
in the literature, we are using d to denote either the struc¬ 
ture map of a differential module, and a non-commutative 
indeterminate). 

Then, if Or is a unit in 21, one can further associate to 
L the quotient 2t(d)/21(d)L ~ 21’’. The differential struc¬ 
ture inherited from 2t(d) makes it a differential module with 
structure map X G X' + CX, where C is the compan¬ 

ion matrix defined above; in other words, this is the module 
(2t(d)/21(d)L, d_c), with the previous notation. 


Scalar extension. Let 21 and IB be two rings equipped 
with derivations and let i/p : 21 —>■ 23 be a ring homomorphism 
commuting with derivation. From a given differential system 
Y' = Ay with coefficients in 21, one can build a differential 
system over 23 by applying p-. it is Y' = p{A)Y, where p{A) 
is the matrix obtained from A by applying p entry-wise. 

This operation admits an analogue at the level of dif¬ 
ferential modules: to a differential module (M, d) over 21, 
we attach the differential module {M<s,d<s) over 23 where 
M<s — 23 M and 9® : M® —>■ M® is defined by: 

V/ € 23, V® G M, d<s{f ® x) = f' ® X f ® d{x). 

It is easily seen that if (M, d) is a differential module asso¬ 
ciated to the system Y' = AY then (M®,9®) is that asso¬ 
ciated to the system Y' = p(A)Y. 


The p-curvature. Let k be any field of characteristic p. We 
assume here that 21 is the field k{x) — consisting of rational 
functions over k — equipped with the standard derivation. 




The p-curvature of a differential module (M, d) over k{x) is 
defined as the mapping : M ^ M. It follows from the 
Leibniz rule (4) and the fact that the p-th derivative of any 
element of k{x) vanishes that the p-curvature is fc(a:)-linear. 

This definition extends to differential systems as follows: 
the p-curvature of the system Y' = AY is the A:(a;)-linear 
map : Ma —>■ Ma where {MA,dA) is the corresponding 
differential module. One can check that the matrix of (in 
the canonical basis of Ma) is the p-th term of the recursive 
sequence (At) defined in (1). 

Considering again a differential operator L and the as¬ 
sociated differential module {'Ql{d)/'Qi{d)L,d-c), for the as¬ 
sociated companion matrix C, we obtain the usual recur¬ 
rence Ai = C and Ai+i = A) + C ■ Ai. The p-curvature of 
2l(9)/2f(9)Z/ will simply be called the p-curvature of L. 

3. SERIES WITH DIVIDED POWERS 

In all this section, we let £ be a ring in which p vanishes. 
We recall the definition of the divided power ring over i, 
and its main properties — mainly, a Cauchy-Lipschitz the¬ 
orem that will allow us to compute solutions of differen¬ 
tial systems. We show how some approaches that are well- 
known for power series solutions carry over without signifi¬ 
cant changes in this context. Most results in this section are 
not new; those from §3.1 and §3.2 are implicitly contained 
in [1, 2], while the theoretical basis of §3.3 is similar to [21]. 

3.1 The ring £[[£1]"^^ 

Let £[[t]]‘^P be the ring of formal series of the form: 

/ = ao -I- ai7i(i) -b a2'y2{t) H--b H- (5) 

where the ai’s are elements of I and each 7 i(t) is a symbol 
which should be thought of as The multiplication on 
£[[t]]‘^P is defined by the rule 7i(i) • 7j(t) = ■ 'li+jit)- 

Remark 3.1. The ring £[[£]]'*'’ is not the PD-envelope in 
the sense of [1, 2] o/£[[i]] with respect to the ideal (t) but its 
completion for the topology defined by the divided powers ide¬ 
als. Taking the completion is essential to have an analogue 
of the Cauchy-Lipschitz Theorem (cf Proposition 3.4). 

Invertible elements of £[[t]]‘^P are easily described: they are 
exactly those for which the “constant” coefficient ao is invert¬ 
ible in £. The ring £[[!]]'**’ is moreover endowed with a deriva¬ 
tion defined by /' = for / = 

It then makes sense to consider differential systems over 
£[[t]]‘^^. A significant difference with power series is the 
existence of an integral operator: it maps / as above to 
/ / = o-iji+iit) and satisfies (f /)' = / for all /. 

Divided power ideals. For all positive integers N, we de¬ 
note by £[[t]]>‘'jv the ideal of consisting of series of the 

form 5I]i>jv The quotient £[[£]]'*'’/£[[f]]t*v i® ^ 

module of rank N and a basis of it is (l, 7 i(t),... , 7 jv-i(i)). 
In particular, for N = 1, the quotient £[[£]]‘^'’/£[[f]]>'i is iso¬ 
morphic to £: in the sequel, we shall denote by /(O) G £ the 
reduction of an element / G £[[£]]‘^® modulo £[[£]]>® . On the 
writing (5), it is nothing but the constant coefficient ao in 
the expansion of /. 

We draw the reader’s attention to the fact that £[[£]]>^ is 

not stable under derivation, so the quotients £[[£l]‘^^/£[[£]]>jv 
do not inherit a derivation. 


Relationship with £[t]. There exists a natural map e : 
£[t] —> £[[£]]‘^^ taking a polynomial atP to aii\-'yi(t). 

The latter sum stops at i = p— 1 because i\ becomes divisible 
by p after that. Clearly, the kernel of e is the principal ideal 
generated by t^. Hence e factors through £[t]/t^ as follows: 

£[t]£[t]/t^ ^ £[[t]f^ (6) 

where pr is the canonical projection taking a polynomial to 
its reduction modulo t^. We observe moreover that the ideal 
t^£[f\ is stable under derivation and, consequently, that the 
quotient ring £[t]/t^ inherits a derivation. Furthermore, the 
two mappings in (6) commute with the derivation. 

3.2 Computations with divided powers 

It turns out that the 7 n(i)’s can all be expressed in terms 
of only few of them, resulting in a more flexible description 
of the ring £[[i]]'*®. To make this precise, we set ti = 7pi(£) 
and first observe that Tf = n\- 7 „j,i it) for all i and n\ this is 
proved by induction on n, using the equalities 

= n\ ■ 7 „j,i(£) • 7 j,i(£) = n! • (£), 

since Lucas’ Theorem shows that 'j = n + 1 (mod p). 

In particular £f = 0 for all i. 

Proposition 3.2. Let n be a positive integer and n = 
X]i=o writing in basis p. Then: 

,no ,ni .n, 

7n(£) = 7no it) ■ 7mp(£) ■ • ■ 7..p« ^ ^ 

no! ni! Us! 

Proof. The first equality is proved by induction on s 
using the fact that if n = a -b &p with 0 < a < p, then 
7a7i)p = 7n, since = 1 (mod p). The second equality 

then follows from the relations t’f' = nd ■ 7„.pi(£). □ 

A corollary of the above proposition is that elements of 
£[[t]]‘^*’ can be alternatively described as infinite sums of 
monomials ang,...,ns ■ £o° ■ ■ ■ ■ t))‘ where the nfs are inte¬ 

gers in the range [0,p) and the coefficient anQ,...,n^ lies in £. 
The product in £[[t]]‘^P is then the usual product of series 
subject to the additional rules if = 0 for all i. 

More precisely, restricting ourselves to some given preci¬ 
sion of the form N = np”, we deduce from the above discus¬ 
sion the following corollary. 

Corollary 3.3. For N = np”, with s G N and n G 
{1,... ,p}, there is a canonical isomorphism of £- algebras: 

£[[t]f^/£[[t]]% cz £[to,...,ts]/it^o,---X-i,t:)- 

For instance, if we take s = 0 and = n in {1,... ,p}, we 
obtain the isomorphism £[[£]]'^^/£[[£]]>^ — £[t]/t^. 

In terms of complexity, the change of bases between left- 
and right-hand sides can both be done in 0~{N) operations 
in £: all the factorials we need can be computed once and for 
all for 0(min(A^,p)) operations; then each monomial con¬ 
version takes 0(s) = 0{log(N)) operations, for a total of 
0(fVlog(fV)) = O'iN). 

The previous corollary is useful in order to devise a mul¬ 
tiplication algorithm for divided powers, since it reduces 
this question to multivariate power series multiplication (ad¬ 
dition takes linear time in both bases). To multiply in 
£[to,..., ts]/ (fp,..., £s_i, £")i one can use a direct algorithm: 
multiply and discard unwanted terms. Using for instance 


Kronecker’s substitution and FFT-based univariate arith¬ 
metic, we find that a multiplication in at precision N 

{i.e. modulo ^[[t]]>^jv) can be performed with 0'(2*°®p 
operations in k. A solution that leads to a cost for 

any e > 0 is in [28], but the former result will be sufficient. 

3.3 The Cauchy-Lipschitz Theorem 

A nice feature of the ring ^[[t]]'^^ — which does not hold 
for £[[t]] notably — is the existence of an analogue of the 
classical Cauchy-Lipschitz theorem. This property will have 
a fundamental importance for the purpose of our paper; see 
for instance [21, Proposition 4.2] for similar considerations. 

Proposition 3.4. Let Y' — AY be a differential system 
of dimension r with coefficients in £[[t]]‘^'’. For all initial 
data V £ F" (considered as a column vector) the following 
Cauchy problem has a unique solution in ^[[t]]'^®’.' 

Y' = A-Y 
T(0) = V. 

Proof. Let us write the expansions of A and Y: 

OO OO 

A = ^ Ai')i{f) and Y = ^ Yi^i{t) 

1=0 i =0 

where the Afs and Yfs have coefficients in £. The Cauchy 
problem translates to Yq = V and Yn+i = Er=o (?) • • 

Yn-i- It is now clear that it has a unique solution. □ 


exponent; then, we have (n — l)p“ < N < np‘. Since we are 
only interested in costs up to logarithmic factors, we may 
assume that we do all operations at precision np" (a better 
analysis would take into account the fact that the precision 
grows quadratically). 

By Corollary 3.3 and the discussion that follows, arith¬ 
metic operations in ■^[[t]]‘^^/^[[t]]>?ips take time 0'(2*°®p ^N). 
This is also the case for differentiation and integration, in 
view of the formulas given in the previous subsection; trun¬ 
cation is free. The total complexity of Algorithm funda- 
mental_solutions is therefore 0~{2^°^p ^Nr‘^) operations 
in i, where r is the dimension of the differential system. If 
N = p^''^\ which is what we need later on, this is 0'{Nr'^). 

The case of differential operators. We now consider 
the case of the differential system associated to a differential 
operator L = ard^-\-- ■ --l-aid-l-ao G £[[t]]‘^'’(9). We will work 
under the following few assumptions: we assume that Or is 
invertible, and that there exists an integer d < p such that 

all fli’s can be written m = ai,o + ai,i 7 i(t) H-h oti^d'yd{t) 

for some coefficients in f; thus, by assumption, Or.o is 
a unit in 1. Our goal is still to compute a basis of solutions 
up to precision N\ the algorithm is a direct adaptation of a 
classical construction to the case of divided powers. 

In all that follows, we let /o,..., f ,— i be the solutions of L 
in f[[t]]‘^'’, such that fi is the unique solution of the Cauchy 
problem (c/ Proposition 3.4): 

i(/0 = 0 ; = forO<j<r (7) 


Of course. Proposition 3.4 extends readily to the case 
where the initial data V is any matrix having r rows. In 
particular, taking V = F (the identity matrix of size r), we 
find that there exists a unique r x r matrix Y with coeffi¬ 
cients in such that T(0) = Ir and Y' = A ■ Y. This 

matrix is often called a fundamental system of solutions. 

Finding solntions nsing Newton iteration. In charac¬ 
teristic zero, it is possible to compute power series solutions 
of a differential system such as Y' = A ■ Y using Newton 
iteration; an algorithm for this is presented on [5, Fig. 1]. 

One can use this algorithm to compute a fundamental 
system of solutions in our context. For this, we first need 
to introduce two notations. Given an element / G ■^[[f]]'^'’ 
written as / = together with an integer m, we 

set [/]"“ = aiji{t). Similarly, if M is a matrix with 

coefficients in ^[[t]]'*’’, we define [M]"* and J M by applying 
the corresponding operations entry-wise. 


Algorithm fundainental_solutions 

Input: a differential system Y' = AY, an integer N 

Output: the fund, system of solutions modulo £[[t]]>?v 

1. y = Ir + t A(0); Z = Ir', m = 2 

2. while m < N/2: 

3. Z = Z +\Z{Ir-YZ)]'^ 

r , “1 2m 

4. y = y- y(/^-(y'-[A]2™-iy)) 

5. m = 2m 

6. return Y 


where Sij is the Kronecker delta. For / = Ej(Lo^f7j (^) 
^[[t]]‘^P, a direct computation shows that the n-th coefficient 
of L{f) is ELoE?=o Assume L{f) = 0. 

Then, extracting the term in fn+r, and using that Oij = 0 
for j > d, we get (r^+r = ^ E[=o E^=o 
Letting m = i-j, we find = EmE_d An.{n)(n+m with 


^ a..o ? 


) CX-i i- 


= E 


, ar,o{i - mf. 

^<r—1 ^ ' 


0< 

0<i — m<.d 


and = n{n — 1) ■ • • (n — (i — m — 1)) is a falling facto¬ 

rial. The expression above for Am is well-defined, since we 
assumed that d < p, and shows that Am is a polynomial of 
degree at most d. 

From this, writing the algorithm is easy. We need two sub¬ 
routines: from_f alling_f actorial(F’), which computes the 
expansion on the monomial basis of a polynomial of the form 
F = X)o<j<n eval]!^, N), which computes the val¬ 

ues of a polynomial F at the N points {0,..., N—1}. The 
former can be done using the divide-and-conquer algorithm 
of [8, Section 3] in time 0~{n); the latter by the algorithm 
of [17, Chapter 10], in time 0''(deg(F) -|- N). The previous 
discussion leads to the algorithm solutions_operator be¬ 
low. In view of the previous discussion, the cost analysis is 
straightforward (at step 2., notice that all required factorials 
can be computed in time 0{d)). The costs reported in the 
pseudo-code indicate the total amount of time spent at the 
corresponding line. 


Correction is proved as in the classical case [5, Lemma 1]. 

Let us take n G {2,..., p} and s G N such that n — 1 is the 
last digit of N written in basis p, and s the corresponding 


Algorithm solutions_operator 

Input: a differential operator I G £[[t]]‘^‘’(cl) of bidegree 

{d, r), with d < p; an integer N 

Output: the solutions fo,..., f ,— i at precision N 






1. for m = —d ,.. ., r — 1: 

9 4 _ 

Cost: 0{d{r + d)) 

3. Am ~ from_falling_f actorial(>lm) 

Cost: 0~{d{r + d)) 

4. Store eval(A„i, A'" — r) 

Cost: 0'{{d + N){r + d)) 

5. for i = 0,..., r — 1: 

6. /i = [0,..., 0,1, 0,..., 0] (ith unit vector of length r) 
Cost: O(r^) 

7. for n = 0,. .., — r — 1: 

8. fi^n^r = ^ Amin) fi,n + m 

Cost: 0{rN{r + d)) 

9. return /o, • • •, fr-i 


Altogether, we obtain the following result, where we use 
the assumption N > d to simplify slightly the cost estimate. 

Lemma 3.5. Suppose thatp < d. Given a positive N > d, 
the classes of fo,. ■ ■, fr-i modulo can be computed 

with at most 0{rN{r + d)) operations in i. 

In particular. Algorithm solutions_operator has a bet¬ 
ter cost than fundainental_solutions when d = 0(r"“^). 


i = k[x]/S (which is thus not necessarily a field). Let a £ £ 
be the class of x. We consider the ring homomorphism: 

ifs : k[x] 

f{x) i-A f{t + a)modt^. 

Regarding the differential structure, we observe that <ps 
commutes with the derivation when £[t]/t''’ is endowed with 
the standard derivation We furthermore deduce from 
the fact that S and /a are coprime that ips extends to a 
homomorphism of differential rings fc[a;][jC] £[t]/t'‘’ that 
we continue to denote by ips- We set ips = l o ips where 
L is the canonical inclusion £[t]/t^ ^ (c/ §3). As 

before, ips commutes with the derivation. Finally, because 
S is separable, we can check that ips is surjective and its 
kernel is the ideal generated by . Hence ps induces an 
isomorphism: 

k[x]/S^^k[x][^]/S^^£[t]/t^. ( 8 ) 

Let Ys be a fundamental system of solutions of the dif¬ 
ferential system Y' = ips (A) ■ Y, i.e. Vs is an r x r ma¬ 
trix with coefficients in such that ys(0) = Ir and 

Ys — Ips (A) ■ Ys- The existence of Ys is guaranteed by 
Proposition 3.4. Moreover, the matrix Ys is invertible be¬ 
cause ys(0) = Ir is. 

Proposition 4.3. Keeping the above notations, we have: 


4. COMPUTING THE P-CURVATURE 

In all this section, we work over a field k of characteristic 
p > 0. We consider a differential system Y' = AY of di¬ 
mension r and denote by Ap the matrix of its p-curvature. 
We write A = ^A, where /a is in k[x\ and A is a matrix 

with polynomial entries. Let d = max(deg /a, deg A), where 
deg A is the maximal degree of the entries of A. We recall 
([13, Prop. 3.2], [9, Lemma Ij) a bound on the size of Ap. 
The bound follows from the recurrence (1), and it is tight. 

Lemma 4.1. The entries of the matrix f^-Ap are all poly¬ 
nomials of degree at most dp. 

The goal of this section is to prove the following theorem. 

Theorem 4.2. There exists an algorithm (presented be¬ 
low) which computes the matrix of the p-curvature of the 
differential system Y' = AY in 0~(pdr‘^) operations in k. 

It is instructive to compare this cost with the size of the out¬ 
put. By Lemma 4.1, the latter is an r x r matrix whose en¬ 
tries are rational functions whose numerator and denomina¬ 
tor have degree ~ pd, so its size is roughly pdr^ elements of k. 
Our result 0~{pdr‘^) is quasi-optimal if we assume that ma¬ 
trix multiplication can be performed in quasi-optimal time. 

4.1 A formula for the p-curvature 

Let Ap denote the matrix of the p-curvature of the dif¬ 
ferential system Y' = AY (in the usual monomial basis). 
The expression of Ap given at the very end of §2 is unfortu¬ 
nately not well-suited for fast computation. The aim of this 
subsection is to give an alternative formula for Ap using the 
framework of divided powers. 

In order to relate k{x) and a ring we pick a sepa¬ 

rable polynomial S £ k\x\ which is coprime with fA and set 


ps{Ap) = -Y^^^ -Ys-^ (9) 

(p) 

where Yg is the matrix obtained from Ys by taking the p-th 
derivative entry-wise. 

Proof. We set Zs = and let (M, d) denote the 

differential module over associated to the differen¬ 

tial system Y' = ips{A)Y. Let yi,...,yr denote the col¬ 
umn vectors of Ys. They are all solutions of the system 
Y' = ips{A)Y, meaning that d{yi) = 0 for all i. Further¬ 
more, if (ei,...,er) is the canonical basis of (^[[t]]'^^)’^, we 
have the matrix relations: *Ys ■ §. = V and e = ^Zs ■ y where 
y (resp. e) is the column vector whose coordinates are the 
vectors yfs (resp. the efs). Applying d to the above rela¬ 
tion, we find 9(e) = *^Zs-y-\-^Zs-d{y) = ^Z's-y and iterating 

this p times, we deduce 9^(e) = ^Z^g'^ ■ y = ^Z^g'^ ■ ‘Ys • e. 
On the other hand, the matrix ips{Ap) of the p-curvature is 
defined by the relation 9^(e) = *‘ips{Ap).e. Therefore we get 
ipsiAp) = Ys ■ Z^g'^. Now differentiating p times the relation 
YsZs = Ir, we find Yg^^Zs -|- Ys • Z^'^ = 0. Combining this 
with the above formula for ips{Ap) concludes the proof. □ 

In our setting, the matrix Ap has coefficients in A:[a:][yj] 
{cf Lemma 4.1), from which we deduce that ips{Ap) has 
actually coefficients in the subring £[t]/t‘‘’ of There¬ 

fore, using Eq. (9), one can compute ips{Ap) knowing only 
Ys modulo the ideal ^[[t]]t 2 p' 

One can actually go further in this direction and establish 
a variant of Eq. (9) giving an expression of ips{Ap) which 
involves only the reduction of Ys modulo ^[[t]]>^. To make 
this precise, we need an extra notation. Given an integer 
i £ [0,p) and a polynomial / £ lltl/i^ (resp. a matrix 
M with coefficients in l\t]/t'^), we write Coeff(/,i) (resp. 
Coeff(M, i)) for the coefficient in T in / (resp. in M). 




Proposition 4.4. Keeping the above notations, we have: 
i^siA,) = -Fs • Fi^^O) • Fg-i 

= Fs-Coeff(^-Fs,p-l)-Fs"^ (10) 
where we have set Ys = Fs mod 

Proof. Differentiating p times the relation Yg = tps{A) ■ 
Ys, we observe that is solution of the same differ¬ 

ential system Y' = 'il)s{A)Y. Hence, thanks to unique¬ 
ness in Cauchy-Lipschitz Theorem, we have the relation 
fJ^^ = Ys ■ fJ^^ (0). The hrst part of the Proposition follows 
by plugging this in Eq. (9) and reducing the result modulo 
To establish the second part, it is now enough to 
notice that the relation Yg = ips{A) ■ Yg implies: 

fJ^’)(0) = (A - = -Coeff(A. Fs, p-1) 

the minus sign coming from (p — 1)! = —1 (mod p). □ 

Remark 4.5. IFe can rephrase Proposition 4-4 o,s fol¬ 
lows: letting yi,... ,yr denote the column vectors of Yg and 
iji G {£[t]/PY be the reduction of yi, the p-curvature of A 
modulo P is the linear endomorphism of {l[t]/PY whose 
matrix in the basis {yi,... ,yr) is Coeff(A • Fs, p—1). It is 
worth remarking that the latter matrix has coefficients in the 
subring £ o/f[t]/t^. 

Remembering Eq. (8), we conclude that Proposition 4.4 
allows us to compute the image of the p-curvature Ap mod¬ 
ulo S^. The strategy of our algorithm now becomes clear: 
we first compute Ap modulo for various polynomials S 
and, when we have collected enough congruences, we put 
them together to reconstruct Ap. The first step is detailed 
in §4.2 just below and the second step is the subject of §4.3. 

4.2 Local calculations 

In all this subsection, we fix a separable polynomial S G 
k[x] and denote by m its degree. Our goal is to design an 
algorithm for computing the matrix Ap modulo . After 
Proposition 4.4, the main remaining algorithmic issue is the 
effective computation of the isomorphism ipg and its inverse. 

Applying ipg and its inverse. We remark that ipg factors 
as follows: 

k[x\/S^ —>■ k[x,t]/{S,{t — xY) —>■ k[x,t]/{S,P) 

X t t + a. 

Applying the right-hand mapping, or its inverse, amounts 
to doing a polynomial shift in degree p with coefficients in 
k[x\/S. Using the divide-and-conquer algorithm of [16], this 
can be done in 0~{p) arithmetic operations in k[x]/S, which 
is 0~{pm) operations in k. Thus, we are left with the left- 
hand factor, say ipg. Applying it is straightforward and can 
be achieved in 0~{pm) operations in k. It then only remains 
to explain how one can apply efficiently . 

We start by determining the image of x by call it 

y = ifig~^{x)-, we may identify it with its canonical preim¬ 
age in k[x\, which has degree less than pm. Write y = 
^o<i<p (i(x’’)x'‘, with every (i in k[x\ of degree less than 
m (so that has degree less than pm). Its image 

through g}*s is Eo<i<pCi(f^)F, which is Eo<i<p Ci(*^)F, 
since = t^ in k[x, t]/{S, {t — x)^). 

Since ipsiy) = x, we deduce that Co(a:^) = x mod S and 
Ciix^) ~ 0 mod S for i = 1,... ,p — 1. The first equality 


implies that x^ generates k[x]/S, so the fact that has 
degree less than m implies that is the unique polynomial 
with this degree constraint such that (o(x^) = x mod S. The 
other equalities then imply that = 0 for i = l,...,p—1. 

In order to compute ((o, we first compute v = x^ mod 
S, using 0'(m log(p)) operations in k. Then, we have to 
find the unique polynomial ((o of degree less than m such 
that Co(i^) = X mod S. In general, one can compute hr 
0(m“) operations in k by solving a linear system. In the 
common case where m < p, there exists a better solution. 
Indeed, denote by tr : k[x]/S —>■ k the fc-linear trace form 
and write U = tr(ir“) and = tT{xK), for i = 0,..., m — 

1. Then formulas such as those in [25] allow us to recover 
(o from t = {to,..., tm-i) and t' = {t'o,..., tY-i) in time 
0~{m). These formulas require that m < p and that S' 
be invertible modulo S, which is ensured by our assumption 
that S is separable. To compute t and t', we can use Shoup’s 
power projection algorithm [29], which takes 0(771^“"*"^^^^) 
operations in k. 

Once ((o is known, to apply the mapping to an el¬ 

ement g{x,t), we proceed coefficient-wise in t. Write g — 
Eo<i<p with all gi of degree less than m. Then 

Ps~^{9) = Eo<i<p ( 5 i(Co) mod T) ( 2 ;J’)F where T is the 
polynomial obtained by raising all coefficients of S to the 
power p, so that S{x)^ — T{x^). 

Computing T takes 0(m log(p)) operations in A:; then, 
computing each term gi{(o) mod T can be done using the 
Brent-Kung modular composition algorithm for 
operations in fc; the total is Finally, the eval¬ 

uation at x^ and the summation needed to obtain gi*g~^{g) 
do not involve any arithmetic operations. 

Remark 4.6. In the case where S = x"* — c (where c £ k 
and p does not divide m), there actually exists a quite simple 
explicit formula for g}'s~^: it takes t to x and x to c‘*x^" 
where n and q are integers satisfying the Bezout’s relation 
pn -I- qm = 1. Using this, one can compute ‘P%~^{g) in 
0~{pm) operations in k in this special case. 

Conclusion. Let us call phiS and phiS_inverse the two 
subroutines described above for computing ipg and its in¬ 
verse respectively. Proposition 4.4 leads to the following 
algorithm for computing the p-curvature modulo . 


Algorithm local_p_curvature 

Input: a polynomial S and a matrix As G Mr{k[x]/S'^) 
Output: the p-curvature of the system Y' = As F 

1. As.^ = phiS(As) 

Cost: 0~{pr^m) operations in k (with m = deg S) 

2. compute a fund, system of solutions Fs G Mr{£[t]/P) 
of the system Y' = As.^F at precision p. 

Cost: 0~{pr‘^) op. in £ using fundamental_solutions 
Remark: Here £ = k[x\/S 

3. Ap,i = Fs • Coeff(Ays, p-l) ■ Y^^ 
at precision 0{t^) 

Cost: 0~{pr‘^) operations in £ 

4. Ap — phiS_inverse(Ap_£) 

Cost: 0~{pr^m‘^) operations in k in general 

0'(pr^m*'“"''^^'^^) operations in fc if m < p 

5. return Ap. 




To conclude with, it is worth remarking that implementing 
the algorithm local_p_curvature can be done using usual 
power series arithmetic: indeed, we only need to perform 
computations in the quotient ^[[f]]‘^^/^[[f]]>p which is iso¬ 
morphic to by Corollary 3.3. Furthermore, we note 

that if we are using the algorithm fmidainental_solutions 
at line 2, then can be computed by performing an ex¬ 
tra loop in fundamental_solutions; indeed the matrix Z 
we obtain this way is exactly 

4.3 Gluing 

We recall that we have started with a differential system 
Y' = AY (with A = r^A) and that our goal is to compute 
the matrix Ap of its p-curvature. Lemma 4.1 gives bounds 
on the size of the entries of Ap. We need another lemma, 
which ensures that we can find enough small “evaluation 
points” (lying in a finite extension of k). Let Fp denote the 
prime subfield of k. 

Lemma 4.7. Given a positive integer D and a nonzero 
polynomial f G k[x], there exist pairwise coprime polynomi¬ 
als Si,Sn € Fp[a:] with n < D such that: 

• EILi deg S'i > D 

• for all i, the polynomial Si is coprime with f and has 
degree at most 1 -f logp(D -I- deg /). 

Proof. Let m be the smallest integer such that p™' > D+ 
deg/. Clearly m< l+log,j(D+deg/) < l-|-logp(D-|-deg/). 
Let Fpm be an extension of Fp of degree m and K be the 
compositum of k and Fpm. Let Si,..., St be the mini¬ 
mal polynomials over Fp (without repetition) of all elements 
in Fpm c K which are not a root of /. We then have 
deg Si < m for all i and Ei=i deg Si > p"* — deg/ > D. 
It remains now to define n as the smallest integer such that 
EiLidegSi > D. Minimality implies EEi^degSi < D 
and thus n < D. Therefore Si,...,S„ satisfy all the re¬ 
quirements of the lemma. □ 

The above proof yields a concrete algorithm for producing 
a sequence Si,..., Sn satisfying the properties of Lemma 
4.7: we run over elements in Fpm and, for each new element, 
append its minimal polynomial over Fp to the sequence (Si) 
unless it is not coprime with /. We continue this process 
until the condition Er=i deg Si > D holds. Keeping in mind 
the logarithmic bound on m, we find that the complexity of 
this algorithm is at most 0~{D + deg /) operations in k. 
Let us call generate_points the resulting routine: it takes 
as input the parameters / and D and return an admissible 
sequence Si,..., S„. 

We are now ready to present our algorithm for computing 
the p-curvature: 

Algorithm p_curvature 

Input: a matrix A written as A = ■ A 

Output: the p-curvature of the differential system Y' = AY 

1. Si, .. . , Sn = generate_points(/yi, d -|- 1) 

Cost: 0~{d) operations in k 

Remark: we have n = 0{d) and deg Si = O(logd), Vi 

2. for i = 1,..., n: 

Ai^p = local_p_curvature(Si, A mod Sf) 

Cost: 0~{pdr‘^) operations in k 

3. compute B G Mr{k\x\) with entries of degree < pd 

such that B = ■ Bi (mod Sf) for all i 


Cost: 0~{pdr^) operations in k 

4. return -E ■ B 

JA 


In view of the previous discussion and Lemma 4.1, the cor¬ 
rectness and the cost analysis of the algorithm p_curvature 
are both straightforward. Hence, Theorem 4.2 is proved. 

We conclude this subsection with three remarks. First, 
when applying Chinese Remainder Theorem (CRT) on line 
3 of Algorithm p_curvature, we notice that all moduli Sf 
are polynomials in x^. This allows the following optimiza¬ 
tion. Writing ■ Bi = (mod Sf (x)) and 

denoting by Cj the unique solution of degree at most d to 
the congruence system: 

Bj{x) = Bi^j{x) (mod Ti(a;)) where ri(x^) = Sf (*) 

we have B = Bj{x)x^. This basically allows us to re¬ 

place one CRT with polynomials of degree dp by p CRT with 
polynomials of degree d. We save this way the polynomial 
factors in log(p) in the complexity. 

Second, instead of working with n polynomials Si, one 
may alternatively choose a unique polynomial S of the form 
S = X"* — a where m > d is an integer not divisible by p 
and a £ k are such that S and /a are coprime. This avoids 
the use of Chinese Remainder Theorem and the resulting 
complexity stays in 0~(pdr‘^) provided that we use Remark 
4.6 in order to compute the inverse of ps- 

Third, we observe that the algorithm p_curvature is very 
easily parallelizable. Indeed, each iteration of the main loop 
(on line 2) is completely independent from the others. Thus, 
they all can be performed in parallel. Moreover, according 
to the first remark (just above), the application of the Chi¬ 
nese Remainder Theorem (on line 3) splits into pr^ smaller 
independent problems and can therefore be efficiently par¬ 
allelized as well. 

4.4 The case of differential operators 

To conclude with, we would like to discuss the case of a 
differential operator L = Ord'" Or-id' ^ -f • • • -I- aid -\- ao 
with Oi G k[x] for all i, of maximal degree d. 

Recall that the p-curvature of L is that of the differen¬ 
tial module {%{&)/%{d)L,d-c), where C is the compan¬ 
ion matrix associated to L as in (3). Applying directly 
the formulas in Proposition 4.4 requires the knowledge of 
the solutions of the system Y' = —CY. It is in fact eas¬ 
ier to compute solutions for the system X' = ^CX, since 
we saw that these solutions are the vectors of the form 
*(?/, y',..., y^' ^^), where p is a solution of L. This is how¬ 
ever harmless: the p-curvatures Ap and Bp of the respective 
systems Y' = —CY and X' = ^CX (which are so-called ad¬ 
joint) satisfy Ap = —‘Bp. Thus, we can use the formulas 
given above to compute ps{Bp), and deduce ps{Ap) for a 
negligible cost. Equivalently, one may notice that the funda¬ 
mental matrices of solutions of our two systems are transpose 
of one another, up to sign. 

Moreover, instead of using the second formula of Propo¬ 
sition 4.4 to compute the local p-curvatures, we recommend 
using the first one, which is (ps(Bp) = —Xs ■ Xg^^(O) • Xg ‘ 
where Xs is a fundamental system of solutions of X' = ‘CX 
and Xs denotes its reduction in Mr{l[t]/t^). If /o,..., f ,— i 
are solutions of the system (7), the (i,/)-th entry of Xs 
is just /j'\ Hence the matrices Xs and Xg^\o) can be 
obtained from the knowledge of the image of ffs modulo 
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Figure 1: Average running time on random inputs of various sizes 


just by reorganizing coefficients (and possibly mul¬ 
tiplying by some factorials depending on the representation 
of elements of we are using). 

As for the fi’s, they can be computed by the algorithm so- 
lutions_operator (provided its assumptions are satisfied). 
We need finally to compute since As(0) is the iden¬ 

tity matrix, this can be done either using Newton iterator, 
a divide-and-conquer approach or a combination of both, 
which computes the inverse of Xs at a small precision, and 
uses divide-and-conquer techniques for higher ones (the lat¬ 
ter being the most efficient in practice). All these remarks 
do speed up the execution of our algorithms when d is not 
too large compared to r. 

Last but not least, we notice that, in the case of differ¬ 
ential operators, the matrix Ap is easily deduced from its 
first column. Indeed, writing Ap = {aij)o<i,j<r and letting 
Cj = -f • • • -I- aijd + aoj G k{x){d) be the dif¬ 

ferential operator obtained from the j-column of Ap, it is 
easily checked that Cj+i is the remainder in the Euclidean 
division of dcj by L. Comparing orders, we further find 
Cj+i = dcj — where lc(cj) is the leading coefficient 

of Cj. This remark is interesting because it permits to save 
memory: indeed, instead of storing all local p-curvatures 
Ap,i, we can just store their first column. Doing this, we 
can reconstruct the first column of Ap using the Chinese 
Remainder Theorem (c/ §4.3) and then compute the whole 
matrix Ap using the recurrence. 

5. IMPLEMENTATION AND TIMINGS 


readable format is about 1GB (after bzip2 compression, it 
decreases to about 300MB). 
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We implemented our algorithms in Magma in the case of 
differential operators; the source code is available at https: //githfilfl 
Figure 1 gives running times for random operators of degrees 
{d,r) in k[x]{d) and compares them with running times of 
(a fraction free version of) Katz’s algorithm which consists 
in computing the recursive sequence (Ai) until i = p. In 
each cell, the first line (resp. the second line) corresponds to 
the running time obtained with our algorithm (resp. Katz’s 
algorithm); a dash indicates that the corresponding running 
time exceeded one hour. Our benchmarks rather well reflect 
the predicted dependence with respect to p: quasi-linear for 
our algorithm and quadratic for Katz’s algorithm. 

Larger examples (than those presented in Fig. 1) are also 
reachable: for instance, we computed the first column of the 
p-curvature of a “small” multiple of the operator con¬ 
sidered in [10, Appendix B.3] modulo the prime 27449. This 
operator has bidegree (d,r) = (108,28). The computation 
took about 19 hours and the size of the output in human- 
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